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Abstract 

We  examine  an  M/M/1  retrial  queue  with  an  unreliable  server  whose  arrival,  service,  failure, 
repair,  and  retrial  rates  are  all  modulated  by  an  exogenous  random  environment.  Provided  are 
conditions  for  stability,  the  (approximate)  orbit  size  distribution,  and  mean  queueing  perfor¬ 
mance  measures  which  are  obtained  via  matrix-analytic  methods.  Additionally,  we  consider  the 
problem  of  choosing  arrival  and  service  rates  for  each  environment  state  with  the  objective  of 
minimizing  the  steady  state  mean  time  spent  in  orbit  by  an  arbitrary  customer,  subject  to  cost 
and  revenue  constraints.  Two  numerical  examples  illustrate  the  main  results. 
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1  Introduction 


This  paper  examines  an  M/M/1  retrial  queue  with  an  unreliable  server  whose  arrival,  service, 
failure,  repair,  and  retrial  rates  are  all  modulated  by  an  external  environment.  For  this  system, 
arriving  customers  who  find  the  server  busy  or  failed,  or  customers  whose  service  is  interrupted 
by  a  server  failure,  join  a  retrial  queue  (or  orbit)  from  which  they  persistently  attempt  to  gain  (or 
regain)  access  to  the  server  at  random  intervals.  Retrial  customers  can  gain  access  to  the  server 
only  when  it  is  found  operational  and  idle,  and  they  repeat  service  until  their  service  requirement 
has  been  satisfied.  The  server  experiences  both  active  and  idle  failures;  the  former  correspond  to 
failures  that  occur  when  the  server  is  processing  a  customer,  and  the  latter  occur  when  the  server 
is  idle.  Once  a  failure  occurs,  the  server  immediately  commences  a  repair  cycle  whose  duration  is 
stochastic.  The  server  cannot  fail  when  it  is  under  repair.  While  the  arrival,  service,  failure,  repair, 
and  retrial  times  are  assumed  to  be  mutually  independent,  they  are  all  modulated  by  a  common 
random  environment.  Like  many  queueing  models  of  this  type,  we  assume  the  environment  is  an 
ergodic  continuous-time  Markov  chain  (CTMC)  on  a  finite  state  space.  We  analyze  this  system 
using  classical  matrix-analytic  methods  (MAM)  developed  originally  by  Neuts  [35],  focusing  on  the 
stability  analysis  and  (numerical)  computation  of  the  steady  state  distribution  of  the  orbit  size 
from  which  we  obtain  steady  state  performance  measures.  Additionally,  we  consider  the  problem 
of  choosing  the  modulated  arrival  and  service  rates  that  minimize  the  mean  time  spent  in  orbit  by 
an  arbitrary  customer  who  arrives  in  steady  state. 

Of  particular  relevance  to  the  present  paper  are  single  server  retrial  models  with  unreliable 
servers  (i.e.,  those  with  servers  that  experience  failures  at  random  intervals).  Three  seminal  papers 
related  to  this  topic  include  [3,  6,  26].  For  retrial  systems  with  an  unreliable  server  and  no  queue 
for  primary  arrivals,  customers  who  arrive  to  find  a  busy  or  failed  server  are  (generally)  routed  to 
the  orbit  in  order  that  they  may  retry  their  service  later  (cf.  [4,  5,  6,  32,  42]).  Other  models  include 
both  a  primary  queue  and  a  retrial  queue  so  that  arrivals  who  are  initially  blocked  from  service 
wait  in  the  primary  queue  (cf.  [7,  17,  18,  31,  37,  39]).  All  of  the  aforementioned  models  assume 
that  the  queueing  system  operates  in  a  static  environment  that  does  not  influence  the  system  in 
any  way. 

The  literature  related  to  retrial  queues  operating  in  random  environments  is  now  beginning  to 
emerge.  Klimenok  [24]  studied  a  BMAP/SM/1  system  whose  operating  mechanism  is  governed 
by  a  random  environment.  In  that  article,  the  limiting  distribution  at  embedded  and  arbitrary 
instants,  as  well  as  the  main  steady  state  performance  measures,  were  examined.  Roszik  and 
Sztrik  [36]  used  the  Modelling,  Specification  and  Evaluation  Language  (MOSEL)  to  analyze  a 
finite-source  retrial  queue  in  a  random  environment  when  all  random  variables  are  assumed  to  be 
exponentially  distributed.  Other  important  models  include  the  BM  AP  /  PH  /I  and  BMAP/PH/N 
systems  analyzed  by  Kim  et  al.  [22,  23]  which  encompass  a  very  broad  class  of  queueing  systems 
with  randomly  varying  rates,  including  the  M/M/1  queue  analyzed  in  [33,  34,  35].  To  analyze  these 
complex  systems,  the  authors  show  that  the  system  state  process  can  be  viewed  as  an  asymptotically 
quasi- Toeplitz  Markov  chain.  Using  results  from  [25],  they  determine  the  stability  condition  and 
devise  an  algorithm  for  computing  the  steady  state  distribution.  Other  recent  contributions  include 
an  examination  of  the  finite-source  MAP/PH/N  retrial  system  with  negative  arrivals  operating  in 
a  random  environment  due  to  Wu  et  al.  [43].  All  of  the  systems  described  here  exhibit  complex 
arrival  and  service  processes;  however,  these  models  do  not  explicitly  consider  the  interplay  between 
a  fully-modulated  system  and  the  impact  of  an  unreliable  server. 

The  model  we  consider  here,  namely  the  unreliable  M/M/1  retrial  queue  in  a  random  envi¬ 
ronment,  could  be  analyzed  as  a  special  case  of  B  MAP /PH /I  retrial  queue  of  [22]  but  for  the 
failure  mechanism  of  the  server.  We  compromise  some  model  complexity  in  the  arrival  and  service 
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processes  with  the  intent  of  (1)  explicitly  accounting  for  an  unreliable  server,  and  (2)  considering 
optimization  of  the  arrival  and  server  rates  for  designing  stable,  efficient  systems  that  meet  quality- 
of-service  guarantees.  Like  most  complex  retrial  queueing  models,  ours  exhibits  the  level-dependent 
quasi-birth-and-death  (LDQBD)  structure  (see  [14,  20,  21,  28]).  We  examine  the  stability  condi¬ 
tion  of  the  system  using  classical  techniques,  namely  Lyapunov  functions.  Furthermore,  we  employ 
matrix-analytic  methods,  via  the  algorithms  of  Bright  and  Taylor  [14,  15],  to  obtain  the  (approx¬ 
imate)  steady  state  distribution  and  important  performance  measures.  Finally,  we  consider  the 
problem  of  choosing  environment-dependent  arrival  and  service  rates  that  minimize  the  mean  wait¬ 
ing  time  in  orbit,  subject  to  a  limited  budget  and  a  minimum  threshold  for  the  revenue  generated 
by  the  system. 

The  remainder  of  the  paper  is  organized  as  follows.  Section  2  provides  a  formal  description  of 
the  retrial  queueing  model  and  the  modulating  environment  and  shows  that  the  queueing  system 
exhibits  the  LDQBD  structure.  Section  3  briefly  discusses  the  form  of  the  limiting  distribution  and 
examines  a  sufficient  stability  condition.  Section  4  reviews  an  algorithm  to  compute  the  steady 
state  distribution  and  mean  performance  measures,  and  provides  a  numerical  illustration.  Finally, 
in  Section  5,  we  consider  the  problem  of  choosing  arrival  and  service  rates  that  minimize  the  steady 
state  mean  time  spent  in  orbit. 

2  Model  Description 

Consider  a  single-server  M/M/1  retrial  queue  operating  in  a  random  environment  whose  server 
is  subject  to  both  active  and  idle  failures.  That  is,  the  server  experiences  failure  whether  it  is  idle 
or  busy,  but  cannot  fail  if  it  is  under  repair.  If  a  primary  arrival  finds  the  server  operational  (i.e., 
not  failed)  and  idle,  it  seizes  the  server  and  immediately  begins  its  service  cycle.  Barring  a  server 
failure  during  this  service  cycle,  the  customer  completes  service  and  departs  the  system.  However, 
should  the  server  fail  during  the  service  cycle,  the  customer  is  immediately  removed  from  the  server 
and  sent  to  an  infinite-capacity  retrial  queue  (or  orbit)  to  retry  service  later.  On  the  other  hand, 
an  arriving  customer  who  hnds  the  server  busy  or  failed  is  routed  directly  to  the  orbit  since  there  is 
no  queue  for  primary  arrivals;  however,  these  customers  are  not  lost.  Retrial  customers  persistently 
attempt  to  regain  access  to  the  server  at  random  intervals,  and  each  behaves  independently  of  the 
primary  arrivals  and  all  other  retrial  customers.  However,  a  retrial  customer  can  only  gain  (or 
regain)  access  to  the  server  if  it  is  found  up  and  idle  at  the  time  of  a  retrial  attempt.  The  service 
discipline  is  preemptive-repeat  (i.e.,  an  interrupted  customer’s  service  cycle  is  repeated  following 
a  successful  retrial  attempt).  All  customers  are  assumed  to  persist  until  they  gain  access  to  the 
server  and  complete  their  service. 

Our  model  is  distinguished  from  other  unreliable  retrial  queueing  models  in  that  its  arrival, 
service,  failure,  repair,  and  retrial  rates  all  vary  randomly  over  time  in  the  spirit  of  the  M/M/1  queue 
in  a  random  environment  studied  by  Neuts  [34].  Specifically,  the  arrival,  service,  failure,  repair,  and 
retrial  rates  are  all  modulated  by  an  external  process  {Z{t)  :  t  >  0}  -  an  irreducible,  continuous-time 
Markov  chain  (CTMC)  with  finite  state  space  S  =  {1, . . .  ,m},  infinitesimal  generator  Q  =  [g^^], 
i,j  G  S,  and  invariant  probability  vector  tt  =  (7ri,7r2, . . .  ,TTm)  that  uniquely  solves  ttQ  =  0  and 
TT  e  =  1  where  0  is  the  zero  vector  of  dimension  m  and  e  is  a  column  vector  of  ones.  When  the 
environment  is  in  state  j  G  S,  primary  customers  arrive  according  to  a  Poisson  process  with  rate 
Xj ,  and  the  service  time  is  exponentially  distributed  with  mean  1/ fij.  When  the  server  is  not  failed 
and  is  either  idle  or  busy,  server  failures  occur  according  to  a  Poisson  process  with  rate  ^j.  Repair 
of  the  server  is  initiated  immediately  following  a  failure,  and  the  duration  of  the  repair  time  (or 
down  period)  is  exponentially  distributed  with  mean  l/aj. 
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For  any  m-dimensional  (row)  vector  x  =  (xi,  X2,  ■  ■  ■ ,  Xm),  denote  its  transpose  by  x',  and 
let  the  diagonal  matrix  of  the  elements  of  x  be  A(a;)  =  diag(xi,  X2,  •  •  • ,  Xm).  Next,  define  the 
m-dimensional  vectors  A  =  (Ai,  A2,  •  •  • ,  A^),  M  fJ-2,  ■  ■  ■ ,  fJ-m),  ^  =  (6, 6,  •  •  •  >  and  a  = 

{ai,a2,  ■  ■  ■  ,am)-  Each  retrial  customer  attempts  to  gain  access  to  the  server  independently  of  all 
other  customers  (primary  or  retrial),  at  exponentially-distributed  time  intervals  with  mean  1/dj 
when  Z{t)  =  j.  Therefore,  if  Z{t)  =  j  and  there  are  i  customers  in  the  orbit,  the  total  retrial  rate 
is  r{i,j)  =  iOj.  Denote  the  vector  of  retrial  rates  by  0  =  (01,02,  •  •  •  ,0m)-  The  row  vectors  A,  /x, 
Q,  and  6  are  all  strictly  positive.  In  the  usual  way,  we  assume  that  the  arrival,  service,  failure, 
repair,  and  retrial  processes  are  mutually  independent;  however,  each  is  modulated  by  the  random 
environment,  {Z{t)  :  t  >  0}. 

For  each  t  >0,  let  R{t)  be  the  orbit  size,  Z{t)  be  the  state  of  the  random  environment,  and  X(t) 
be  the  status  of  the  server  so  that  X{t)  =  0  if  the  server  is  failed,  X{t)  =  1  if  the  server  is  up  and  idle, 
and  X{t)  =  2  if  the  server  is  up  and  busy  at  time  t.  The  state  of  the  queueing  system  is  described 
by  the  continuous-time  stochastic  process,  {R,Z,X)  =  {{R(t) ,  Z (t) ,  X (t))  :  t  >  0},  with  state 
space  E  =  {{i,j,  k)  :  i  >  0,  j  G  S,  k  £  {0, 1,  2}}.  Note  that  E  contains  one  denumerable  dimension 
(the  orbit  size)  and  two  hnite  dimensions  (the  environment  state  and  server’s  status).  Because  all 
inter-arrival,  service,  inter-failure,  repair,  and  inter-retrial  times  are  exponentially  distributed,  it  is 
easy  to  see  that  {R,  Z,X)  is  a  continuous-time  Markov  chain  (CTMC)  on  E.  Proposition  1  asserts 
that  this  CTMC  has  a  well-structured,  block  diagonal  infinitesimal  generator  matrix. 


Proposition  1  The  process  {R,  Z,X)  with  state  spaee  E  is  a  level- dependent  quasi-birth- and- death 
(LDQBD)  process  with  block  diagonal  infinitesimal  generator  matrix 
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where  Ci  =  Q  —  A(A  -|-  ^)  —  iA{6),  i  >  0,  Di  =  Q  —  A(A  -|-  cc),  and  D2  =  Q  —  A(A  -|-  p  -|-  ^). 


The  LDQBD  process  (cf.  [14,  21,  28])  is  a  natural  extension  of  the  QBD  process  wherein  some  or 
all  of  the  matrices  comprising  the  fth  level  are  explicitly  dependent  on  the  level  i.  For  our  purposes 
here,  the  limiting  distribution  of  the  LDQBD  process  is  needed  to  compute  the  steady  state  queueing 
performance  measures  and  to  formulate  and  solve  optimization  problems  to  determine  optimal  (or 
near  optimal)  operating  rates  for  each  of  the  m  distinct  environment  states.  Assuming  its  existence, 
define  the  limiting  distribution  of  {R,  Z,X)  as  the  row  vector  p  =  (Po,Pi,P2,  •  •  •)  where  for  f  >  0, 
Pj  is  a  3m-dimensional  row  vector  of  limiting  probabilities  restricted  to  level  i.  If  the  process 
is  ergodic,  p  is  the  unique  positive  solution  to  pQ*  =  0  and  pe  =  1.  A  closed- form,  necessary 
and  sufficient  condition  for  the  positive  recurrence  of  general  LDQBD  processes  is  not  available; 
however,  it  is  known  that  if  p  exists,  then  it  has  the  matrix-geometric  property. 
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Let  us  assume  for  momentarily  that  {R{t) ,  Z (t) ,  X (t))  — )■  {R,Z,X)  in  distribution  as  t  — )■  oo. 
In  such  a  case,  p  =  where  p{i,j,k)  =  F{R  =  i,  Z  =  j,  X  =  k),  (i,j,k)  £  E.  The 

marginal  distribution  of  the  steady  state  orbit  size  is  given  by 

m  2 

P(^  =  i)  =  ^  =p.e,  z>0, 

j=l  k=0 

and  likewise,  the  steady  state  status  of  the  server  has  probability  mass  function  (p.m.f.) 

OO  m 

F{X  =  k)  =  '^'^p{i,j,k),  fcG  {0,1,2}. 

i=Q  j  =  l 

In  Section  3,  we  provide  the  stability  condition  of  the  queueing  system  using  Lyapunov  functions 
to  establish  a  sufficient  condition  for  the  ergodicity  of  {R,  Z,X). 


3  Stability  Analysis 


In  this  section  we  discuss  necessary  and  sufficient  conditions  for  positive  recurrence  of  the 
process  {R,  Z,X).  The  following  theorem  can  be  stated  using  results  in  Bright  and  Taylor  [14]. 


Theorem  1  The  LDQBD  process  {R,Z,X)  with  infinitesimal  generator  matrix  Q*  is  positive  re¬ 
current  if  and  only  if  there  exists  a  strictly  positive  solution  to  the  system  of  equations 


Po  (^^0  +  Rq  ©i)  —  0) 


subject  to  the  normalization  condition 
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In  such  a  case,  the  3m-order  row  vector  is  given  by 

2—1 


Pi = Po  n  ^ 


n=0 


(2) 

(3) 

(4) 


In  equations  (3)  and  (4),  when  i  =  0,  the  empty  product  results  in  the  identity  matrix  I.  It  is  well 
known  (cf.  [14,  28])  that  the  sequence  {Ri  :  i  >  0}  is  the  minimal  non- negative  solution  of  the  set 
of  equations 

A  +  i?i  Tj+i  +  Ri  (i?j_|_i0j_|_2)  =  0,  i  >  0  (5) 

which  must  be  determined  numerically. 

In  general,  it  is  difficult  to  assert  positive  recurrence  using  the  conditions  of  Theorem  1.  For  our 
retrial  queueing  system.  Theorem  2  establishes  a  sufficient  condition  for  stability  using  Lyapunov 
functions  and  a  classical  result  due  to  Tweedie  [41],  which  is  stated  here  as  Lemma  1. 


Lemma  1  A  continuous-time  Markov  chain  with  generator  matrix  Q*  =  [(Z^j./],  x,x'  £  E  is  regular 
and  ergodic  if  there  exists  a  function  v  :  E  ^  which  is  bounded  below,  a  finite  set  H  C  E,  and 
some  e  >  0  such  that 


and 


d{x)  =  q*,,,/[v{x')  —  u(x)]  <  oo,  for  all  x  £  E, 

x'eE\{x} 

d{x)  <  — e,  for  all  x  £  E\H . 
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Equipped  with  Lemma  1,  we  are  now  prepared  to  state  a  sufficient  condition  for  the  ergodicity 
of  the  continuous-time  process  {R,  Z,X).  Theorem  2  represents  a  somewhat  strong  condition. 

Theorem  2  The  process  {R,Z,X)  is  ergodic  if 

A' <  [A(q:  +  ^)]“^  A(q;)/i',  (8) 

where  the  inequality  holds  componentwise. 


Proof.  Since  all  the  inter-event  times  are  assumed  to  be  exponential,  and  the  environment  is 
a  CTMC,  we  consider  the  following  Lyapunov  function.  For  i  >  0,  j  £  S,  and  some  real  numbers 
a,  6  >  0,  define 

b  +  ai,  if  fc  =  0, 


v{i,j,  k) 


1  +  ai,  if  fc  =  1, 


_az,  if  fc  =  2. 


Substituting  the  function  v  in  (6),  and  using  the  transition  rates  in  Q*  of  (1),  we  obtain 


d(i,  j,  0)  =  a \j  —  OLj  b, 

d{i,  j,  1)  =  \j  +f,jb  +  idj{l  -  a), 

d{i,  j,  2)  =  a  Xj  -  pLj  +  f,j{b  -  I  +  a). 

It  is  not  difficult  to  see  that  d{i,j,  k)  <  oo  for  all  {i,j,  k),  so  we  need  only  to  determine  if  the  drift 
functions  satisfy  (7).  Clearly,  for  all  j  £  S,  there  exists  a  positive  integer  N  such  that  d{i,j,  1)  <  0 
for  each  f  >  A  if  a  >  1.  Therefore,  inequality  (7)  is  verified  if  there  exist  positive  constants  a  and 
b  satisfying  the  set  of  linear  inequalities. 


1  —  a  <  0, 

(9) 

aXj  —  Ujb  <  0, 

(10) 

aXj  —  fij  +  f,j{b  —  1  +  a)  <  0. 

(11) 

By  inequality  (10),  we  must  have  b  >  aXjfaj.  Using  this  fact  in  inequality  (11),  we  can  eliminate 
b  to  obtain 

^3^3  +  +  ^3^3 

But  by  (9),  we  must  have  a  >  1;  therefore,  the  set  of  linear  inequalities  (9)-(ll)  has  a  solution  if 
and  only  if  there  is  a  positive  number  a  such  that  a  >  1  and 


^  ^  +  hj) 

Xj{aj  +f,j)+  ’ 

or  equivalently,  if 

^3<{  j  =  l,2,...,m.  (12) 

Rewriting  this  expression  in  vector/matrix  form,  we  obtain  <  A{a.  A{a)fj,' ,  where  the 

inequality  holds  componentwise.  ■ 

The  right-hand  side  of  (12)  can  be  viewed  as  the  effective  service  rate  when  the  environment 
is  in  state  j  as  the  quantity  ajj{aj  +  Cj)  is  the  effective  proportion  of  time  the  server  is  not  failed 
in  environment  state  j.  Likewise,  Xj  is  the  effective  arrival  rate  of  customers  in  state  j.  Under 


6 


the  strong  condition  that  the  effective  arrival  rate  is  less  than  the  effective  service  rate  for  all 
environment  states,  stability  is  to  be  expected;  however,  condition  (12)  need  not  hold  for  every 
j  G  S.  To  see  this,  note  that  the  average  effective  arrival  rate  is  given  by  A  =  tvX' ,  as  the  arrival 
process  mirrors  that  of  the  standard  M/M/1  queue  in  a  random  environment  (see  Neuts  [34]). 
Similarly,  using  a  Markov  reward  argument,  it  can  be  shown  that  the  average  effective  service  rate 
is  given  by  /r  =  7rA(Q;  +  Standard  results  for  stability  conditions  of  single-server 

retrial  queues  (see  [8,  19])  show  that  A  <  /r  is  necessary  for  stability  of  the  system.  Hence,  the 
system  can  be  stable  only  if 

ttA'  <  7rA(Q: -|- ^)~^  A(q;)/xA  (13) 

But  it  is  possible  that  Xj  >  aj  fJ-j/ioij  +  for  some  j  G  S  while  (13)  still  holds.  Therefore,  we  see 
that  while  (12)  implies  (13),  it  is  not  necessary  for  stability  of  the  system.  For  convenience,  and 
for  use  in  Sections  4  and  5,  let  us  define  the  overall  traffic  intensity  by 


ttA^ 

^  7rA(Q:  +  4)“^A(Q:)/r' 


(14) 


Remarks:  Note  that  if  ^  =  0,  the  server  never  fails,  and  the  stability  condition  reduces  to 
<  TT/r'.  This  is  precisely  the  stability  condition  for  the  standard  M/M/1  queue  in  a  random 
environment  analyzed  by  Neuts  [34].  Moreover,  if  Xj  =  A,  gij  =  ^j  =  and  aj  =  a  for  all 
j  G  S,  the  (necessary  and  sufficient)  stability  condition  is  A  <  a[i/{a  +  The  same  result  has 
been  derived  for  other  single-server,  exponential  retrial  models  with  an  unreliable  server  (cf.  Falin 
[17]),  or  it  can  be  derived  from  general  service  time  models  by  assuming  exponential  service  times 
(cf.  Sherman  et  al.  [38,  39],  Kulkarni  and  Choi  [26],  and  others). 


4  Evaluating  Performance  Measures 

The  most  widely  used  algorithms  for  approximating  the  limiting  distribution  of  a  LDQBD 
process  with  an  infinite  number  of  levels  or  phases  are  due  to  Bright  and  Taylor  [14,  15].  Because 
these  algorithms  play  a  central  role  in  solving  the  optimization  problem  of  Section  5,  we  next 
summarize  their  essential  elements.  These  techniques  extend  the  well-known  logarithmic  reduction 
algorithm  of  Latouche  and  Ramaswami  [27]  QBD  processes  to  the  level-dependent  case. 

The  generator  matrix  of  (1)  possesses  two  nice  properties  that  facilitate  relatively  easy  imple¬ 
mentation  of  the  algorithms.  First,  the  number  of  phase  states  in  each  level  is  fixed  at  3m,  and 
second,  the  matrix  A  is  independent  of  the  level  i.  Essentially,  the  main  algorithm  of  [14]  truncates 
the  infinite  series  of  (3)  at  some  level  K  and  then  re-normalizes  to  compute  an  approximate  sub¬ 
vector  of  the  form  Pi{K)  =  Pq{K)  On^o  ^  >  1,  where  Pq{K)  satisfies  (2)  and  the  normalization 


condition,  Po(-^)  0^=0'^"  e  =  1.  The  subvectors,  {pi{K)  :  i  >  0},  represent  an  invariant 

measure  for  the  limiting  distribution  of  all  states  at  or  below  level  K;  therefore,  p^  <  Pi{K)  compo¬ 
nentwise  for  any  K  >  0,  and  Pi{K)  — >■  p^  componentwise  as  iF  — )>  oo.  For  a  given  truncation  point 
K,  Bright  and  Taylor  [14]  examine  the  discrete-time  Markov  chain  embedded  at  the  jump  epochs 
of  the  process  to  obtain  the  family  of  matrices  {Ri  :  *  >  0}  using  a  recursive  scheme.  Lemma  2  is 
a  direct  consequence  of  (1)  and  Lemma  1  of  [14]. 

Lemma  2  If  {R,Z,X)  is  positive  recurrent,  the  matrix  Ri  is  given  hy 


T*-l 


r-i 


i£— 1— n 
2+2^-^  ’ 


i  >  0 


(15) 


f=0 


n=0 
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where  for  i  >  1,  and  Df  are  3m  x  3m  matrices  recursively  defined  by 


=  A(-F*+i 

)-\ 

II 

CD 

T 

T' 

ir\ 

ur 

''  ^ i+2‘+^  ^ i+Z-2^  ^i+2^+'^'-^i+2‘ 

—  F)^ 

T  F)^  F)^  ifi 

The  infinite  series  of  equation  (15)  can  be  truncated  using  a  simple  scheme  (see  Algorithms  2  and 
3  of  [14]).  By  rearranging  the  terms  in  (5),  the  family  of  matrices,  {Ri  :  i  >  0},  are  recursively 
computed  by 

Ri  =  ^  (— Tj+i  —  Ri+i  0j+2)  ^  (16) 

(assuming  the  inverse  exists)  so  that  (15)  need  not  be  computed  repeatedly.  For  all  of  the  numerical 
results  that  follow,  Algorithms  1-3  of  [14]  are  used  to  select  the  integer  K  and  compute  pfiK). 

Assuming  p  <  1,  the  approximate  steady  state  performance  measures  of  the  queueing  system 
{R,Z,X)  may  be  obtained  using  the  approximation  of  the  steady  state  distribution  given  by  the 
vector  p  =  k)],  where  k)  =  P(i?  =  i,  Z  =  j,  X  =  k),  {i,j,  k)  G  E.  More  specifically,  the 

steady  state  orbit  size  distribution  is 


m  2 

P(^  =  i)  =  p(i,j,A:)  =p.e,  z>0.  (17) 

j=l  k=0 

Likewise,  the  steady  state  distribution  of  the  server’s  status  is 

oo  m 

7fe  =  P(X  =  fe)  =  A:),  A:  =  0,1,2.  (18) 

j=0  j=l 

(Of  course,  the  distribution  of  Z  is  the  invariant  vector  tt,  which  is  independent  of  {R,X).) 

From  (17)  and  (18),  we  can  obtain  the  steady  state  delay  and  congestion  parameters  in  the 
usual  way  using  Little’s  Law.  Specifically,  using  the  (approximate)  distribution  p,  the  steady  state 
mean  orbit  size  is 

OO 

nR)  =  Y.^{p^e).  (19) 

i=l 

Let  N  denote  the  steady  state  number  of  customers  in  the  system  (in  the  retrial  queue  and  in 
service).  The  mean  of  N  is  easily  computed  by  noting  that  N  =  R  ii  X  =  0orX  =  1,  and 
A^  =  i?  +  lifX  =  2.  Therefore,  by  conditioning  on  X,  we  see  that 

E(iV)  =  E{R)  (70  +  71)  +  E(i?  +  1)72  =  E{R)  +  72.  (20) 

We  note  that  the  mean  number  in  system  is  not  E{R)  +  p  because  there  are  periods  during  which 
the  orbit  is  not  empty,  but  the  server  is  idle.  Next,  let  W  be  the  sojourn  time  (time  in  service 
and  in  orbit)  of  an  arbitrary  customer  who  arrives  in  steady  state.  It  is  well  known  that,  for  an 
ordinary  (non- modulated) ,  single-server  retrial  queue  with  Poisson  arrivals  and  exponential  inter¬ 
retrial  times,  the  mean  sojourn  time  is  the  mean  number  in  system  divided  by  the  arrival  rate. 
Analogously,  we  can  apply  Little’s  Law  to  obtain 

E(W)  =  (7rA')“^E(7V), 
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(21) 


where  (ttA^)  ^  is  the  average  effective  arrival  rate  of  customers  to  the  system.  In  a  similar  manner, 
the  steady  state  mean  time  spent  in  the  orbit  is  given  by 


^{Wr)  =  (7rA')“^E(i?) 


(22) 


Finally,  to  compare  various  environment  states,  we  define  the  traffic  intensity  in  state  j  by 

-1 


Pj  ~ 


“i  + 


(23) 


Next,  we  provide  a  numerical  example  to  illustrate  the  steady  state  orbit  size  distribution, 
among  other  steady  state  performance  measures. 


Example:  Suppose  the  environment  has  state  space  S  =  {1,  2,  3, 4,  5,  6,  7}  and  infinitesimal  genera¬ 
tor  matrix 


-7.6 

2.0 

3.0 

1.0 

1.0 

0.1 

0.5 

0.5 

-8.0 

2.5 

1.0 

2.0 

0.8 

1.2 

0.3 

1.5 

-5.8 

1.0 

1.0 

1.0 

1.0 

2.0 

3.0 

5.0 

-11.2 

0.1 

1.0 

0.1 

0.8 

2.5 

2.0 

1.1 

-7.9 

0.7 

0.8 

1.5 

1.0 

1.6 

1.2 

0.5 

-6.3 

0.5 

2.0 

2.5 

2.0 

1.0 

1.8 

0.9 

-10.2 

Table  1  reveals  that  the  system  is  critically  loaded  when  the  environment  is  in  state  3  (pa  = 
0.96),  and  it  is  overloaded  in  environment  state  2  (p2  >  !)•  The  environment  spends  a  relatively 
short  amount  of  time  in  the  overloaded  condition  (less  than  20%)  and  about  30%  of  the  time  in 
a  critically  loaded  condition.  However,  the  overall  traffic  intensity  is  only  p  =  0.6895.  On  first 
glance,  this  result  seems  counterintuitive  until  we  notice  the  relatively  small  arrival  and  failure 
rates  of  state  2.  In  environment  state  3,  the  service  rate  is  comparatively  high,  so  even  though  the 
system  is  critically  loaded  in  this  state,  the  system  will  experience  recovery  periods  when  occupying 
the  less  detrimental  states. 


Table  1:  Summary  of  parameter  values  and  traffic  intensities. 


Environment  (j) 

A. 

Pj 

% 

Pj 

P 

1 

3.0 

7.0 

0.5 

2.0 

1.0 

0.1021 

0.5357 

0.6895 

2 

1.0 

3.0 

1.1 

0.5 

11.0 

0.1917 

1.0667 

3 

3.0 

12.5 

1.5 

0.5 

2.0 

0.3086 

0.9600 

4 

2.0 

12.5 

4.0 

2.0 

2.0 

0.0848 

0.4800 

5 

2.0 

4.5 

1.0 

6.0 

5.0 

0.1256 

0.5185 

6 

0.5 

2.0 

0.7 

3.0 

1.0 

0.1130 

0.3083 

7 

0.5 

4.0 

1.5 

0.5 

0.5 

0.0740 

0.5000 

The  values  of  the  approximated  performance  measures  are  summarized  in  Table  2.  It  is  note¬ 
worthy  that  the  limiting  probability  that  the  server  is  busy  (P(X  =  2))  is  less  than  the  traffic 
intensity  p.  This  is  attributed  to  the  fact  that  there  are  periods  in  which  the  retrial  queue  is  not 
empty,  but  the  server  remains  idle  waiting  for  either  a  primary  customer  arrival  or  the  next  retrial 
attempt. 
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Table  2:  Performance  measures  for  the  numerical  example. 


E{R) 

E{N) 

E{Wr) 

E(VF) 

F{X  =  0) 

P(X  =  I) 

cT 

II 

5.9904 

6.2963 

3.0902 

3.2480 

0.4685 

0.2256 

0.3059 

For  the  sake  of  completeness,  Figure  1  graphs  the  (approximate)  probability  distribution  of  R, 
namely  F{R  =  i)  =  Pj  e  for  i  =  0, 1, , . . . ,  75.  Note  the  geometric  rate  of  decay  exhibited  by  the 
steady  state  distribution,  which  is  expected  because  p  is  a  matrix  geometric  distribution. 
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Figure  1:  Approximate  steady  state  orbit  size  distribution. 


In  Section  5  we  use  the  approximated  steady  state  orbit  size  distribution  and  queueing  perfor¬ 
mance  measures  to  determine  optimal  arrival  and  service  rates  to  minimize  the  mean  time  spent 
in  orbit. 

5  Optimizing  Arrival  and  Service  Rates 

In  this  section,  we  consider  the  problem  of  choosing  the  arrival  and  service  rates,  as  a  function 
of  the  environment  state,  to  minimize  the  mean  time  spent  in  orbit  by  an  arbitrary  customer  in 
steady  state.  This  problem  can  be  viewed  as  a  (static)  design  problem  in  which  the  system  designer 
chooses  Aj  G  [A,  A]  and  fij  G  [m,  jz]  for  each  j  G  S  where  0  <  A  <  A  <  oo  and  0  <  fj.  <  JI  <  oo. 
The  rate  setting  is  done  only  once  so  that  whenever  Z(t)  =  j,  a  controller  limits  the  arrival  rate  to 
the  optimal  Xj  value  and  tunes  the  service  rate  to  the  optimal  fij  value,  j  =  1,  2, . . . ,  m.  For  every 
unit  of  arrival  rate,  the  system  gains  a  reward  (revenue)  rj,  and  for  every  unit  of  service  rate,  the 
system  incurs  a  cost  Cj.  Let  r  =  (ri,  r2, . . . ,  r^,)  and  c  =  (ci,  C2, . . . ,  Cm)  be  the  revenue  and  cost 
vectors,  respectively.  The  total  revenue  generation  rate  (over  all  environment  states)  must  meet  or 
exceed  a  minimum  threshold  value  R  {0  <  R  <  oo),  while  the  total  cost  rate  is  subject  to  an  upper 
limit  B  {0  <  B  <  oo).  It  is  reasonable  to  assume  that  only  admission  and  service  rates  can  be 
chosen  at  our  discretion  as  the  failure  and  repair  rates  are  dictated  by  the  inherent  limitations  of 
the  equipment,  etc.  Similarly,  the  retrial  rates  may  be  dictated  by  customer  behavior  or  attitudes 
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and  are  assumed  to  be  outside  of  the  controller’s  purview. 

The  objective  is  to  minimize  the  mean  time  spent  in  orbit  by  an  arbitrary  customer  in  steady 
state  given  by 

OO  OO 

=  (24) 

i=l  i=l 

Here,  we  include  the  implicit  dependence  on  /x  through  the  vectors  p^,i>0,  which  are  approximated 
by  the  algorithms  summarized  in  Section  4.  To  compute  'd(A,  /x),  we  truncate  the  infinite  series  in 
(24)  at  the  nth  term  (n  G  N)  if  |S'n+i  —  Sn\  <  e  where 

n 

Sn  =  e), 

i=l 

and  e  >  0  is  a  convergence  threshold.  With  these  preliminaries  and  notation,  the  nonlinear  pro¬ 
gramming  formulation  is  as  follows; 

min  'd(A,  p) 


P<1, 

(25a) 

r\'  >R, 

(25b) 

cp'  <  B, 

(25c) 

Aj  G  [a,  a]  , 

.  ,m, 

(25d) 

Pj  G  [p,  p]  , 

i  =  1,  •  ■ 

. ,  m. 

(25e) 

Constraint  (25a)  enforces  the  necessary  stability  condition  discussed  in  Section  3  while  (25b)  and 
(25c)  are  the  revenue  and  cost  constraints.  Constraints  (25d)  and  (25e)  are  box  constraints  that 
ensure  realistic  rate  settings.  For  a  candidate  solution  (A,  p),  the  left-hand  side  of  (25a)  is  computed 
directly  via  (14). 

Though  easily  stated,  the  optimization  problem  (25)  is  not  easily  solved  for  at  least  the  fol¬ 
lowing  reasons.  First,  the  objective  fnnction  is  expensive  as  it  requires  an  approximation  of 
p  =  iPQ,Pi,P2,  ■  ■  ■)  for  each  candidate  solution  (A, /x).  Second,  the  feasible  region  is  not  closed 
(due  to  the  strict  inequality  constraint  (25a)),  and  in  general,  derivative  information  about  the  ob¬ 
jective  function  is  not  available.  These  complications  motivate  our  use  of  derivative-free,  adaptive 
search  algorithms  with  proven  convergence  properties.  We  next  discuss  a  class  of  these  algorithms. 

5.1  Solving  the  Rate-Setting  Problem 

Because  the  steady  state  vector  p  =  (pg ,pi,p2,  ■  ■  ■)  is  only  available  numerically  (via  the  matrix- 
analytic  methods  described  in  Section  4),  the  objective  function  (24)  is  computationally  expensive. 
Therefore,  to  solve  problem  (25),  we  employ  generalized  pattern  search  (GPS),  and  specifically, 
derivative- free,  mesh-adaptive  search  (MADS)  techniques.  These  techniques  do  not  require  a  closed- 
form  objective  function,  as  long  as  the  objective  function  can  be  evaluated  numerically  at  points 
inside  the  feasible  region.  Detailed  descriptions  of  these  algorithms  are  snmmarized  in  [2,  12,  13]. 

GPS  is  a  derivative- free  optimization  technique  for  unconstrained  problems  originally  intro¬ 
duced  by  Torczon  [40]  who  proved  convergence  of  a  subsequence  of  iterates  to  a  first-order  station¬ 
ary  point.  It  has  known  convergence  properties  for  a  variety  of  problem  classes,  even  when  the 
objective  function  is  nonsmooth  (see  Audet  and  Dennis  [11]).  GPS  methods  iteratively  search  a  set 
of  points  around  the  cnrrent  iterate  for  one  that  improves  the  objective  function  value.  Gonsider 
the  general  nonlinear  minimization  problem, 

min??(®),  (26) 


11 


where  G  :  £  <  Ax  <  u],  d  :  M”  — )•  M  and  A  G  ig  a  rational  matrix.  Moreover,  £ 

and  u  are  the  lower  and  upper  bounds  of  the  constraints  where  £,  w  G  {M™  H  {±oo}}  and  i  <  u. 
GPS  algorithms  generate  a  sequence  of  iterates  {x^}  in  M"'  with  nonincreasing  objective  function 
values.  Each  iteration  is  divided  into  an  optional  search  step  and  a  local  poll  step.  Both  the 
search  and  poll  steps  evaluate  points  on  a  mesh  in  order  to  find  a  mesh  point  that  improves  the 
objective  function  value.  The  mesh  is  constructed  as  a  lattice  of  points  in  M"’,  based  on  a  finite  set 
of  directions  D  that  form  a  positive  spanning  set  and  a  mesh  size  parameter,  (A^  >  0),  that 
controls  the  fineness  of  the  mesh.  In  this  case,  a  positive  spanning  set  refers  to  a  set  of  vectors 
such  that  any  vector  in  the  space  can  be  represented  by  a  nonnegative  linear  combination  of  the 
vectors  in  the  set.  By  definition,  nonnegative  linear  combinations  of  the  elements  of  the  set  D  span 
M"’.  The  directions  that  form  D  can  be  arbitrarily  chosen  provided  that,  for  each  direction  dj  G  D, 
j  =  1,  2, . . . ,  |Z)|,  dj  =  Gzj,  where  G  G  is  a  nonsingular  matrix  and  Zj  G  Z”  is  an  integer 

vector.  At  iteration  k,  the  mesh  is  centered  around  the  current  iterate  x^  G  M”  and  its  fineness  is 
parameterized  through  the  mesh  size  parameter  A^.  The  mesh  can  then  be  represented  as 

Mk  =  ^Xk  +  AfcD  z  :  z  G  z[|_  ^  (27) 

where  Z+  is  the  set  of  nonnegative  integers.  Note  that  in  (27),  the  columns  of  the  matrix  D  form 
the  set  D. 

In  the  search  step,  GPS  can  evaluate  any  finite  set  of  mesh  points,  and  a  number  of  strategies 
exist  for  generating  trial  points,  including  random  search,  genetic  algorithms,  Latin  hypercube 
search,  or  orthogonal  arrays.  If  the  search  step  fails  to  provide  an  improved  mesh  point,  the  poll 
step  is  invoked.  The  poll  step  is  more  rigidly  defined  and  evaluates  the  neighboring  mesh  points  of 
the  current  iterate.  The  use  of  positive  spanning  directions  in  the  construction  of  these  neighboring 
points  provides  the  theoretical  basis  for  the  convergence  of  GPS.  The  poll  set  at  iteration  k  can  be 
expressed  as  {xk  +  A^d  :  d  G  D^},  where  C  L)  is  also  a  matrix  whose  columns  positively  span 
M"'.  The  poll  set  is  therefore  composed  of  mesh  points  neighboring  the  current  iterate  Xk  in  the 
directions  of  the  columns  of  D^,  a  multiple  A^  away  from  the  current  iterate. 

If  the  search  and  poll  step  both  fail,  the  incumbent  solution  is  said  to  be  a  mesh  local  optimizer 
and  the  mesh  is  then  refined  by  setting  the  mesh  size  parameter 

Afc+i  =  V’“''=Afc,  (28) 

where  >  1  is  rational  and  Wk  G  {w~  ,w~  +  1, . . . ,  —1}  for  some  w~ .  An  incumbent  point  Xk  is 
replaced  by  Xk+\  only  if  ^^{xk+i)  <  '&{xk),  and  in  such  case,  Xk+i  is  termed  an  improved  mesh  point. 
If  an  improved  mesh  point  is  found  in  either  step,  then  the  mesh  is  either  retained  or  coarsened 
by  increasing  the  mesh  size  parameter  according  to  equation  (28)  for  some  Wk  G  {0, 1, . . . ,  u)^}.  It 
follows  that  for  any  k  >  0  there  exists  an  integer  rk  £  such  that  A^  = 

The  convergence  analysis  of  pattern  search  is  well-established  in  [10,  40]  and  requires  a  few 
assumptions.  First,  all  iterates  produced  by  GPS  must  lie  in  a  compact  set  [16].  This  very  common 
assumption  holds  as  long  as  {tc  G  :  d^x)  <  d(a;o)}  is  compact.  Second,  if  the  matrix  G  =  I 
(as  is  usually  the  case),  then  the  constraint  matrix  A  must  be  rational.  The  final  necessary 
assumption  is  that  d(xo)  <  oo  for  xq  G  M”.  Torczon  [40]  proved  that,  under  these  assumptions,  the 
mesh  size  parameter  satisfies  liminffc_,.oo  A^.  =  0,  which  leads  to  the  main  directional  convergence 
result  (see  [11]).  However,  Audet  [9]  proved  convergence  to  a  Clarke  first-order  stationary  point 
in  the  one-dimensional  case  for  unconstrained  problems.  Of  particular  relevance  to  our  work  here, 
GPS  was  extended  in  [29,  30]  to  problems  with  bound  and  linear  constraints,  respectively.  To 
handle  these  constraints  while  maintaining  convergence  properties,  infeasible  points  are  discarded 
without  being  evaluated,  and  search  directions  are  chosen  so  as  to  conform  to  the  geometry  of 
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the  nearby  constraint  boundaries.  The  NOMADm  optimization  software  by  Abramson  [1],  written 
in  the  Matlab®  computing  environment,  was  used  to  implement  the  pattern  search  procedure 
described  here.  NOMADm  is  specifically  designed  to  numerically  solve  nonlinear  and  mixed  variable 
optimization  problems  via  an  implementation  of  the  class  of  mesh-adaptive  direct  search  (MADS) 
algorithms.  GPS  is  a  subclass  of  MADS  [12],  in  which  poll  directions  are  restricted  to  a  uniformly 
bounded  finite  set. 


5.2  Optimization  Illustration 

In  this  subsection,  we  formulate  and  solve  an  instance  of  problem  (25)  using  generalized  pattern 
search  (namely  MADS).  The  vectors  a,  6,  r,  and  c  are  specihed  and  fixed  (i.e.,  they  are  simply 
treated  as  parameters  in  the  model).  The  decision  variables  are  the  arrival  and  service  rates  con¬ 
tained  in  A  and  p,  respectively.  Initial  feasible  vectors  Aq  and  /Xq  were  specified  for  each  of  three 
independent  replications  to  ensure  that  the  algorithm  produced  consistent  solutions.  The  aim  is 
to  choose  the  arrival  and  service  rates  that  minimize  the  mean  time  customers  spend  in  orbit  in 
steady  state. 


Example'.  Suppose  the  retrial  system’s  environment  has  state  space  S  =  {1,2, 3, 4,  5}  and  inhnites- 
imal  generator  matrix 


Q  = 


-23 

6 

6 

4 

8 


6 

-19 

2 

1 

1 


7 

1 

-18 

3 
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9 

4 

9 

-11 

1 


1 

8 

1 

3 

-13 


The  invariant  probability  vector  of  Q  is  tt  =  (0.1946,  0.1071,  0.1697,  0.3530,  0.1754)  while  the 
revenue  threshold  value  is  i?  =  6,  and  the  upper  limit  on  the  budget  is  i?  =  20.  The  box  constraints 
(25d)  and  (25e)  for  this  example  are  \j  G  [1,  5]  and  p.j  G  [0,  4],  j  G  S.  The  other  input  parameters 
are  as  follows:  ^  =  (0.5, 1.1, 1.5, 4.0, 1.0),  cx  =  (2.0,  0.5, 8.5, 4.5,  6.0),  6  =  (1.0, 4.0,  2.0,  8.0,  5.0), 
r  =  (1.0, 1.0,  0.5,  2.0, 0.75),  and  c  =  (2.0, 1.5, 0.5,  2.0,  0.5). 

Table  3  summarizes  the  best  obtained  solutions  using  three  distinct  initial  feasible  solutions. 
The  objective  function  values  of  runs  1  and  3  are  very  similar,  as  are  their  solutions  with  the 
exception  of  the  first  two  elements  of  /x*.  Run  2  produced  a  solution  that  differs  significantly  from 
the  others  and  yields  a  clearly  inferior  objective  function  value.  The  average  number  of  iterations 
needed  for  convergence  of  the  MADS  algorithm  was  just  over  400  (i.e.,  on  average,  400  approximated 
steady  state  distributions  were  computed). 


Table  3:  MADS  best  obtained  solutions  for  the  example  problem. 


Run  no. 

Initial  solution 

Best  solution  obtained 

^{\*,PL*) 

1 

Ao  =  (I.0,1.0, 1.0, 1.0, 1.O) 

/Xq  =  (2.0,  2.0,  2.0,  2.0,  2.0) 

A*  =  (1.008, 1.287, 1.001, 1.211, 1.046) 
PL*  =  (1.917,3.983,3.965,3.105,3.999) 

30.791 

2 

Ao  =  (1.0,1.0, 1.0, 1.0, 1.0) 
/xo  =  (1.5, 1.5, 1.5, 1.5, 1.5) 

A*  =  (1.102, 1.263, 1.000, 1.171, 1.057) 
PL*  =  (2.389,3.416,3.741,3.476,2.552) 

33.293 

3 

Ao  =  (1.0,1.0, 1.0, 1.0, 1.0) 
/Xq  =  (2.5,  2.5,  2.5,  2.5,  2.5) 

A*  =  (1.083, 1.409, 1.001, 1.123, 1.017) 
PL*  =  (3.403,1.962,3.950,3.139,3.995) 

29.580 
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It  is  worth  noting  that  the  time  and  computational  effort  to  solve  the  5-state  example  greatly 
exceeded  the  time  and  effort  needed  to  solve  a  similar  3-state  model  which  required  just  under  100 
iterations.  It  is  surmised  that  the  increased  computational  effort  stems  from  the  exponential  in¬ 
crease  in  the  effort  needed  to  compute  Pq,Pi,P2,  •  ■  •  at  each  objective  function  evaluation.  Despite 
the  computational  effort,  the  solutions  were  obtained  in  less  than  10  minutes.  These  optimization 
models  can  be  used  to  help  a  system  controller  determine  the  appropriate  set  of  arrival  and  service 
rates  to  choose,  depending  on  the  prevailing  conditions. 
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